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ABSTRACT 

The size distribution of asteroids and Kuiper belt objects in the solar system is difficult to reconcile with a bottom-up formation 
scenario due to the observed scarcity of objects smaller than ~100 km in size. Instead, planetesimals appear to form top-down, with 
large 100 - 1000 km bodies forming from the rapid gravitational collapse of dense clumps of small solid particles. In this paper 
we investigate the conditions under which solid particles can form dense clumps in a protoplanetary disk. We use a hydrodynamic 
code to model the interaction between solid particles and the gas inside a shearing box inside the disk, considering particle sizes 
from sub-millimeter-sized chondrules to meter-sized rocks. We find that particles down to millimeter sizes can form dense particle 
clouds through the run-away convergence of radial drift known as the streaming instability. We make a map of the range of conditions 
(strength of turbulence, particle mass-loading, disk mass, and distance to the star) which are prone to producing dense particle clumps. 
Finally, we estimate the distribution of collision speeds between mm-sized particles. We calculate the rate of sticking collisions and 
obtain a robust upper limit on the particle growth timescale of ~10^ years. This means that mm-sized chondrule aggregates can grow 
on a timescale much smaller than the disk accretion timescale (~10^ - 10^ years). Our results suggest a pathway from the mm-sized 
grains found in primitive meteorites to fully formed asteroids. We speculate that asteroids may form from a positive feedback loop in 
which coagualation leads to particle clumping driven by the streaming instability. This clumping, in turn reduces collision speeds and 
enhances coagulation. Future simulations should model coagulation and the streaming instability together to explore this feedback 
loop further. 
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1. Introduction 


Planetesimals are super-kilometer-sized bodies that are the seeds 
of terrestrial planets as well as the cores of ice and gas giants 


( Safronov|1972[[Qiiang & Youdin|2010} [Johansen et al.|2Q14| ). 

Hence, the production of planetesimals is an important step in 
planet formation. In the solar system, large asteroids and Kuiper 
belt objects are the left-over planetesimals that did not become 
incorporated into a planet. Similar planetesimal belts are known 
to exist around other stars ( |Wyatt|2008] ). 


The size distribution of the solar system planetesimals (i.e. 
large asteroids and Kuiper belt objects) reveals a “knee” at 
around 100 km in size, with a relative scarcity of bodies smaller 
than -100 km. Some authors have argued that this knee can¬ 
not be reproduced by bottom-up coagulation scenarios, unless 
the minimal initial size of a planetesimal was also -100 km. 
Hence, it appears that the size of solids in the protoplanetary disk 
“jumped” from the sub-meter scale to -100 km without passing 
through the intermediate sizes, and that bodies smaller than 100 
km formed later by collisional grinding ( [Morbidelli et al.|2009| ). 
A similar bump is inferred for the trans-Neptunian population 


(Fraser & Kavelaars|2009l|Sheppard & Trujilloj2010 

Shankmanj 

et al. 2013t Fraser et al. 2014|). That said, Weidenschi 

ing (|2010|) 


suggests an alternate scenario in which the “knee” is produced in 
a bottom-up process from an initial population of 0.1 km-sized 
bodies. In this scenario, the bump at -100 km is produced by 
the transition from dispersion-dominated runaway growth to a 
regime dominated by Keplerian shear. 


There are also important theoretical difficulties in building 
planetesimals in a bottom-up process. The first hurdle is that 
small particles in the mm-cm size range do not easily stick to 
form larger objects. Indeed, the largest grains observed in a pro 
toplanetary disk are mm-cm in size ( |Testi et ar]|2014t 


van der 


|Marel et~ari|2013| ). To the extent that meter-sized objects do 
form, they quickly spiral into the central star due to friction with 
gas orbiting the star at sub-Keplerian speeds ( [Weidenschilling 


1977}. Axisymmetr ic pressure bumps may stop p article drift ([Jo¬ 


hansen et al.|2009a| ), but growth rates remain low ( [Johansen et al 
2008| ). For example, bodies at 3 AU can grow to a maximum of 


100 m in 1 Myr ( j Windmark et al. |20 1 2| ) . 


These constraints suggest a different picture of how planetes¬ 
imals form. It begins with the growth of macroscopic particles 
by coagulation, followed by the accumulation of these particles 
by a hydrodynamical process like the streaming instability (see 
below) which increases the local concentration of solids. 


y _ ^solid _ ^solid ^ ^ x 

T V-lV’ 

^total ^gas ^ solid 

where Ssoiid and Sgas are the surface densities of the solid and gas 
components of the disk respectively. The solids also sediment 
toward the mid-plane of the disk, so that the volume density of 
solids vs gas at the midplane can be significantly higher than Z. 
If the density of solids reaches the Roche density. 
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where Q is the Keplerian frequency, then the particle self-gravity 
overwhelms Keplerian shear, leading to gravitational collapse 
( [Goldreich & Ward] 1973 1 ). This allows a collection of small ob¬ 
jects to coalesce directly into fully formed planetesimals, and 
explaining the large birth sizes. 

The main obstacle to achieving high densities is disk turbu¬ 
lence, thought to be caused by the magneto-rotational instability 
(MRI, |B albus & Hawley 1 1991 1 ) . But even inside an MRI-free 
“dead zone” ( |Gammie||1996| ), the Kelvin-Helmholtz instability 
can be very effective at producing turbulence that blocks par¬ 
ticle sedimentation to the mid-plane before they can reach the 
Roche dens ity ( |Weidenschilling| 1 9^ |Youdin & Shu|2002[ |Lee 


et al.||2Q10| ). Particle densities high enough for a local gravita¬ 


tional collapse can be obtained by concentration in long-lived 
gaseous vortices ([Barge & Sommeria|1995[ Meheut et al.|2012|). 


or pressure bumps caused by the MRI (Jo 


lansen et aL||2009a 


2Q11| ). These regions can quickly accumulate meter and decime 


ter sized bodies, and for meter-sized boulders, the concentrations 
can significantly reduce radial drift caused by friction with the 
surrounding gas ( [Johansen et al.|20()6||2Q07| ). 

The streaming instability is another powerful mechanism to 
concentrate particles in localized dense clumps and to prevent 
radial drift ([Youdin & Goodman|20Q5[[Johansen et al.|2007[[Jo-[ 


hansen & Youdin|2007[ Bai & Stone|2010a[ ). It is driven by the 

relative drift between the solid and gas components of the disk. 
The gas component of the disk experiences a radial pressure sup¬ 
port that counters the gravitational force of the star. This leads 
to a speed difference Av between the gas and solid components 
of the disk. A useful way to understand the streaming instabil¬ 
ity is that solids experience a head wind from the gas, and the 
gas is, in turn, pushed forward by solids. A small overdensity of 
solids will have a stronger back-reaction on the gas, and hence, 
a lower radial drift than neighboring regions. The reduced radial 
drift leads to the further accumulation of solids as neighboring 
particles drift into the overdensity. In this way, the streaming 
instability bears some similarity to a traffic jam. 

The behavior of solid particles embedded in a gaseous ac¬ 
cretion disk is primarily determined by the Stokes number, Tf = 
ti Qk, which measures the particle friction time tf in terms of the 
Keplerian frequency Qk- Some authors have raised the concern 
that the streaming instability has mainly been studied for rela¬ 
tively large particles, with Tf ~ 0.1, corre sponding to dm-sized 
particles when applied to the asteroid belt ( [Shi & Chiang|2Q13] ). 
Particles of this size are inconsistent with coagulation experi¬ 
ments which show that particles cannot grow beyond mm sizes 
( Zsom et al.|201Q[ [Guttler et al.|2010[ ). Therefore, the main goal 
of this paper is to study the streaming instability for particle sizes 
down to millimeter. 

Chondrules make up most of the mass in primitive mete¬ 
orites, with a typical size of R ~ 0.3 mm, and large chondrules 
reaching ~ 1 mm (e.g. [ J acquet[20 1 4[ ) . At this size range, chon¬ 
drules couple with the small-scale turbulent eddies in the disk, 
and they concentrate in the high-pressure regions between the 
eddies ( [Cuzzi et^[2001[ ). This alone does not lead to particle 
growth, because a typical collision between two chondrules in 
the disk will result in bouncing ( [Guttler et al.|2010| ). However, 
a collision between a small /i m-sized dust particle and a chon- 
drule does result in sticking. [Ormel et al.[ ([2008J have shown 


that chondrules acquire a porous dust layer around them, which 
absorbs some of the kinetic energy of collisions and allows the 
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chondrules to stick more easily. Under typical disk conditions, 
most of the chondrule and dust mass is in the form of millimeter¬ 
sized aggregates. For low turbulence regions, the chondrule ag¬ 
gregates have typical sizes around ~ 4 mm, and maximum 
sizes as high as R ~ 10 cm ( [Ormel et al.[[2008[ ). Our work es¬ 
tablishes the connection between the streaming instability and 
chondrules and chondrule aggregates. While we test the lim¬ 
its of the streaming instability at both the small-particle and the 
large-particle limits, the bulk of our analysis focuses on the size 
range of chondrules and chondrule aggregates. 

In addition to particle size, the streaming instability also de- 
pends on the magnitude of the radial pressure gradient ([Bai & 


Stone|201Q^, as well as the solid concent ration ([Johansen et alT 


2009b|[Bai & Stone|201Qa[ ). In particular, [Bai & Stone[ ( [2Q10a l 


have shown that the streaming instability is effective in produc¬ 
ing particle clumps for particle size Tf > 0.1, and solid concen¬ 
tration Z > 0.02. 

A number of global processes increase the local solid con- 


centration, including radial drift of 

Darticles from the outer disk 

into the inner disk ( 
bumps or vortices 

Youdin & Shu 

[2002[), large scale pressure 

Johansen et a 

. [2009a[), and lastly, there 


is a region near the ice line where ice-dust aggregates (“dirty 
snowballs”) break up on a timescale comparable to the infall 
timescale, leading to a s ignificant accu mulation of solids by as 
much as a factor of 6.7 ( [Sirono[[2011[ ). The second way to in¬ 
crease Z is to remove gas from the disk, either through late-stage 
photoe vaporation (e.g . Alexander et al.|2006a|b[) or thr ough disk 
winds ( [Bai & Stone [201 3[ [Suzuki & 


nutsuka 


[2009[ ). Gas re¬ 


moval has the double effect of increasing Tf and Z simultane¬ 
ously. Hence, gas removal may be a powerful way to trigger the 
streaming instability in the late stages of the disk. 

This paper is structured as follows. In section we present 
our numerical model, including the simulation setup and initial 
conditions. In sectionj^we present our core result - a map of the 
region in the Z vs Tf phase space that is consistent with particle 
clumping (Fig.[^. In section]^ we review the implications and 
offer new constraints on planetesimal formation models (Fig.[^. 
In section we measure the collision speed for Tf = 0.003 par¬ 
ticles and find that particle growth occurs faster than the disk 
accretion timescale. In sectionwe present our conclusions. 


2. Numerical model 

We use the Pencil Code ( [Youdin & Johansen[[2007[ ) to model a 
small vertical slice of the protoplanetary disk in which the solar 
system planets formed using a stratified model that includes ver¬ 
tical gravity. The canonical model for this disk is the minimum 
mass solar nebula (MMSN, [Hayashi[[l981[ ). In it, the surface 
density of the gas component of the disk follows the power law 


9 / r \-3/2 

2;= 1700gcm-2(—] . (3) 

In addition to the gas, the disk also contains solid particles that 
are initially /im in size and represent about 1% of the disk mass. 
Over time, particle sizes can grow by coagulation, and the mass 
ratio between gas and solids may change if particles migrate 
through the disk, or if the gas becomes depleted. 


2.1. Aerodynamics of solid particles 

The gas component of the disk produces a drag force on the par¬ 
ticles whenever the particles have a non-zero relative speed Vrei 
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with respect to the gas. The form of the gas drag law depends 
on the size of the particle relative to the mean free path A of 
gas particles. Small particles, with size R < 9A/4 experience 
Epstein drag, while larger particles experience Stokes drag. We 
will show in section [2^ that the the smaller particles in our sim¬ 
ulations lie in the Epstein regime. Particles in the Epstein regime 
experience the drag force 


With Eq. this gives the particle Stokes number in terms of the 
local surface density 


Tf = 


p.Rn 


( 8 ) 


Combined with Eq. we obtain the relation between R and Tf 
for the MMSN 


An 2 

-Edrag — ^ ^ ^th t^rel 


(4) 


where p is the gas density, R is the particle size, and Vth is the 
mean thermal speed of gas molecules. 


Vth 


8kr 

( TTyUmp 


Here k is the Boltzmann constant, T is the local temperature, p 
is the mean molecular weight of the gas, and is the proton 
mass ( |Armitage|2010[ ). It is convenient to write Vth in terms of 
the isothermal sound speed Cs, 


Note that p, = 3.5 g cm ~^ is a typical density for a mm-sized 
chondrule ( [Hughes|1980| ), and r = 2.5 AU lies in the inner part 
of the present day asteroid belt. Note also that this equation is 
only valid within the Epstein regime. In section [2^ we show that 
this corresponds to Tf < 1. Our two largest particle sizes, Tf = 3 
and 10, lie in the Stokes regime. These particles are a few meters 
in size. 

2.3. Pressure support 



The friction time 


The gas component in a protoplanetary disk experiences a radial 
pressure support that partially counters the gravitational force 
from the central star. As a result, the gas in the disk orbits at 
sub-Keplerian speed. The difference between Keplerian speed 
Vk and and gas speed can be written as 


^ _ mVrel 

^drag 

measures how long it takes for a particle with mass m to have its 
relative speed changed by order unity as a result of gas drag. Eor 
a particle in the Epstein regime (Eq.[^ with material density p,, 
we can write the friction time as 


Av = Vk - = 77 Vk, 


where the pressure gradient parameter rj is given by 

l/csf^lnP 
^ 2 hk / 5 In r ’ 


( 10 ) 


( 11 ) 


p,R 

tf = - 

PVth 


p.R 


(5) 


where p is the gas density. We commonly express tf in terms of 
the Keplerian frequency Qk to obtain the Stokes number 


where P is the gas pressure ( [Nakagawa et al.|1986| ). In the ab¬ 
sence of gas, the solid component of the disk would orbit at the 
Keplerian speed Vk. In the presence of gas, the solids experience 
a “headwind” from the gas component, with a wind speed of Av. 
This headwind is a source of drag that can lead to radial drift. 
Eor scale-free numerical simulations, it is helpful to normalize 
the headwind speed by the sound speed and write 


p.R 

Tf = i2k tf = i2k 

PCs 

The Stokes number is a scale-free measure of the stopping time. 
Since the disk scale height is ^ the Stokes number can 

be written as 



p.R [n 

Jh ys' 


( 6 ) 


2.2. Size of particles in the Epstein regime 

We are interested in the formation of planetesimals at the mid¬ 
plane of the disk, where the gas density is 


Av Vk 1 / Cs \ 5 In P 
Cs ^ Cs 2\vk/51nr’ 


( 12 ) 


This A is a free parameter that measures the strength of the pres¬ 
sure support and the headwind. The pressure is P = pc^. Equa¬ 
tion [7] gives the gas density at the midplane, which combined 
with Qk = c^/H allows us to write 


Ellk c. 


H 


-kCs 

oc r Cs. 




The isothermal sound speed Cs is determined primarily by the 
local temperature 


H 


(7) 



60 m 




(13) 
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where k is the Boltzmann constant, T is the temperature, 
is the proton mass, and fi is the mean molecular weight. In 
the MMSN the disk is assumed to be optically thin, so that 
its temperature profile is T = 280 K(r/However, 
|Chiang & Youdin| ( |2010| ) proposed a colder disk, with T = 
For the moment, then, we will simply write 
the temperature profile as 


-PI2 




Therefore, the pressure P oc r ^ follows a power law with 
exponent 


d\nP 

d\nr 




Equation can now be written as 



(14) 

Using the parameters from the MMSN and the Chiang-Youdin 

nebula (CY) we get 


/ r 


Ammsn - 0-068 (2 5 , 

(15) 

/ r \2/7 


^ = »“'’(2.5Au) ■ 

(16) 


Given the uncertainty in the disk model, we choose A = 0.05 
as a typical value for r ~ 2.5 AU (present day asteroid belt). 
This is also the value chosen by |Bai & Stone| ( |2010a| ), so that 
using A = 0.05 facilitates comparison with their work. We also 
ran simulations with A = 0.025 to study the effect of a reduced 
pressure support. This can occur as part of a large -scale pressure 
bump, which could form around the snow line (|Kretke & Lin 
[20071 ). 

2.4. Epstein vs Stokes regimes 

Particles with R < 9/1/4 experience Epstein drag, where the 
mean free path A is given by 


2.5. Simulation setup 

In our simulations we use Cs = flk = 1 as a natural scaling unit 
for the problem. The free parameters for the problem are A and 
Tf given by Eqs. [^and| |T^ We use the Pencil Cod^to simulate 
the dynamics of gas and solid particles in a 2-dimensional ver¬ 
tical slice of the protoplanetary disk ( jYoudin & Johansen|2007| ). 
We model a box of dimensions 0.2H x 0.2H in the radial-vertical 
plane centered on the midplane of the disk. The box has periodic 
boundary conditions. Our baseline simulations have a 128 x 128 
grid resolution with 16,384 “super-particles” representing the 
solids. Each super-particle represents a large number of phys¬ 
ical particles. We run a total of 102 baseline simulations as we 
vary the key simulation parameters: 

- We test seventeen particle sizes spaced equally in log scale 
from Tf = 0.001 to Tf = 10. 

- We test both the canonical pressure support A = 0.05, and a 
reduced pressure support A = 0.025. 

- We run each (Tf, A) simulation three times, with different 
initial particle positions. The particle positions are chosen 
from a random distribution, uniform within the box. 

Every simulation begins with a solid fraction of Z = 0.005 
(Eq. 0- First, we allow the solids to sediment to the midplane. In 
most cases, the sedimentation is complete within 50 orbits. For 
the Tf < 0.003 simulations we allowed 250 orbits for the sedi¬ 
mentation phase. Once the particles reach an equilibrum scale 
height, we begin the numerical experiment: We reduce the gas 
density p exponentially, so that it halves every 50 orbits. We keep 
Tf constant, which is equivalent to reducing the particle size R at 
the same rate as p (see Eq. [^. This procedure is also equiva¬ 
lent to increasing the particle density Pp while keeping p and R 
fixed. The objective of the experiment is to determine the value 
Z needed to produce particle clumps for a given Tf. 


2.6. System response time 

It is desirable that the rate of gas removal be slower than the 
response time of the system. Otherwise our procedure will tend 
to over-estimate the Z value needed to produce particle clumps. 
The response time of the system is in the order of the particle 
crossing time Across. which is 


0.2H 


(17) 


where 0.2H is the width of the box and Vj- is the particle radial 
velocity, which is given by 


A = 


ncr 


where cr = 2x 10 cm^ is the collision cross section of an H 2 
molecule, and n is the particle number density, which is 


pmp 


Using Eq.[7] along with H = Cs/Ok and Eq.[^ one can compute 
the mean free path. For the MMSN, the mean free path is T ~ 
43.6 cm at 2.5 AU. That means that particles smaller than R ~ 
98 cm lie in the Epstein regime. From Eq.j^we find that all our 
particles with Tf < 1 lie in the Epstein regime. 

Article number, page 4 of[^ 
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-2CsA 
Tf + t:' 


O.IH^ -K 

Across — . vf ^f 


Cc A 


(18) 


(19) 


With this one can show that 4ross < 50 orbits o Tf -k t^^ < 50;:, 
for A = 0.05. This roughly corresponds to the interval 0.0064 < 
Tf < 157. This means that for most of our simulations the system 
response time is less than 50 orbits, but for Tf < 0.0064 we will 
overestimate the Z needed for particle clumps. For Tf < 0.0064, 
our results should be considered robust upper limits. This is a 
necessary limitation because these simulations also suffer from 
small time steps which force long computation times. 


^ https://code.google.eom/p/pencil-code/ 
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3. Results 

A reduced pressure support (A = 0.025) is associated with 
slower radial drift, lower gas turbulence, and a correspondingly 
higher concentration of solids at the midplane. We can quantita¬ 
tively measure the scale height of particles (ifp) and the degree 
of turbulence using the root-mean-square of the particle posi¬ 
tions (z coordinate) and vertical speed. 



where Zi is the z coordinate of the particle, and is the verti¬ 
cal component of the particle’s velocity. Figureshows Hp and 
Vrms across the entire range of simulations. For particles smaller 
than Tf ~ 0.3, the lower A is associated with both a lower scale 
height and lower vertical motion. Particles larger than Tf ~ 0.3 
are poorly coupled with the gas, and are largely unaffected by A. 
The effect of particle size and concentration is more complex. 

We found it helpful to divide the particles into three differ¬ 
ent size bins, based on their qualitative behavior: (1) Particles 
smaller than Tf ~ 0.003 form a suspension, where particles can¬ 
not decouple from the gas even at high concentration and particle 
clumps do not occur. (2) Between Tf ~ 0.003 and Tf ~ 0.3 is the 
optimal range for particles to engage in the streaming instability. 
(3) Particles larger than Tf ~ 0.3 suffer strong radial drift, and the 
over-densities that do form are quickly destroyed by gas erosion. 


3.1. Suspension regime (Vf < 0.003) 

Our smallest particles have Tf = 10“^ and Tf = 10“^-^^ 0.002. 

These particles do not sediment easily because they are strongly 
coupled to the gas. Therefore, we wait 250 orbits for the solids 
to sediment before we begin the experiment. Figure]^ shows the 
evolution of Hp, including the sedimentation phase. The long 
sedimentation phase ensures that Hp has converged before we 
start removing gas. 

Figurej^shows that, for Tf < 0.003, the particle size has very 
little effect on Hp. Instead, Hp is determined by A. For A = 0.05, 
the particle concentration is also important for Z > 0.04. Some 
authors have proposed that, for well-coupled particles such as 
these, there is a critical Z beyond which turbulent mixing be¬ 
comes ineffective, resulting in a high-density “cusp” at the mid¬ 
plane ( |Sekiya|1998 |Youdin & Shu|2002] ). The critical value of 
Z depends on the disk conditions, and range from Z ~ 0.02 for 
a cold disk, and Z ~ 0.10 for a MMSN dYoudin & Shu|[2002]) . 
To test this idea, we ran 256 x 256 simulations with Tf = 10“^-^. 
The particle density can be written as 


Pp(z) = 




( 22 ) 


where z is the vertical coordinate. Figure shows the midplane 
particle density pp(z = 0) as a function of Z, and a snapshot of 
Pp(z) at Z = 0.08. Although we find an increase in pp(z = 0) 
with Z, it i s not nearly of the magnitude proposed by |Youdin &| 
|Shu| ( |2002| ). We suspect that our simulations have an additional 
source of mixing that interferes with the formation of the Sekiya- 
Youdin-Shu cusp. For r = 1 AU, our 128 x 128 grid cells have a 


physical size around 5,000 - 8,000 km. Both the 128 x 128 and 
256 X 256 runs should be able to resolve the cusp seen in Fig. 
1 of |Youdin & Shu|f2002| ). On the other hand, the discrepancy 
between the 256 x 256 runs and the 128 x 128 runs indicates 
that our results have not converged, and higher resolutions may 
reveal higher densities. 


3.2. Streaming regime (0.003 < Tf < 0.3) 

After Tf ~ 0.003, solid particles begin to decouple from the gas, 
leading to a regime where the radial drift is sufficient to make the 
streaming instability effective. Figures]^ andshow spacetime 
diagrams for a selection of runs across all regimes, from Tf = 
10“^ to Tf = 3. For a standard pressure gradient (A = 0.05), there 
are visible, long-lived filaments for 0.003 < Tf < 0.3. For Tf = 1 
the filaments are short lived. Somewhere between Tf = 0.3 and 
Tf = 1 particle clumps start to become difficult again. 

When A = 0.05, the smallest particles that form visible fil¬ 
aments are have Tf = 0.003. The pile-up process is somewhat 




Fig. 1. Particle scale height H^ (top; see Eq. 1^, and root-mean- 
square of the vertical component of the particle speed Vj-ms (bottom; see 
Eq.[^ at the end of the sedimentation phase. The sound speed is set to 
Cs = 800 m s"k Half of the simulations (red squares) have a standard 
pressure support (A = 0.05) and the other half (blue triangles) have 
a reduced pressure support (A = 0.025). Eor particles smaller than 
Tf ~ 0.3, lower A results in stronger sedimentation (lower Hp). 
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analogous to a traffic jam: An initial overdensity of solids has 
a stronger back-reaction on the gas. This reduced A, lowers the 
radial drift, leading to a pile-up as other solids (with faster radial 
drift) enter the filament and increase the density even further. 
Somewhere between Tf = 0.03 and 0.1, clumping starts to be¬ 
come difficult again, as the smaller solid cross section (fewer 
particles, higher mass per particle) cannot lower A as effectively. 
After around Tf > 1, traffic jams are no longer easily observable. 

3.3. Radial drift regime (zf >0.3) 

The regime with Tf > 0.3 is not a good environment for as¬ 
teroid formation. Although particles in this regime do clump, 
they continue to rapidly drift toward the star, and the clumps that 
form are very transient as they are easily eroded by wind in the 
disk. Figure shows the radial speed for various particle sizes 
and concentrations. Higher concentrations of solids are associ¬ 
ated with a decrease in radial speed, but even at high Z the radial 
speed remains high. In this regime, larger particles are more able 
to resist gas drag and exhibit lower radial drift. However, even at 
Tf ~ 10, the radial drift is as high as |Vr| ~ O.OlCg. At that speed, 
these particles would cross the disk in around a thousand years. 

Als note that the overdensities that do form in this regime 
are very short lived, as they are easily destroyed by erosion from 
the gas component. But in any case, it seems unlikely that this 
regime is at all important for realistic disks. As noted earlier, the 
“particles” in this regime are large. A Tf = 10 particle is actually 
a 5-m eter boulder. The inefficiency of sticking ( [Guttler et al.| 
|2010| ), and the streaming regime, probably present formidable 
barriers for any solids to reach the radial drift regime. 

3.4. Conditions for particle streaming 

To locate particle clumps we look for stable peaks in the solid 
surface density I^soiid- The linear stability analysis of jYoudin"^ 
[Goodman] ( |2005| ) shows that the streaming instability grows only 
in modes which have both a radial and a vertical variation. The 
non-linear simulations of [Johansen & Youdin] ( [200'?] ) also show 
significant vertical variation in the particle density. However, 
simulations which include the vertical gravity on the particles are 


stratified; they have a dense midplane layer and little additional 
variation in the vertical direction. For this reason, the vertical av¬ 
erage captures very well the magnitude of particle overdensities 
in the stratified simulations that we present here. 

Figure illustrates our strategy. We are interested in iden¬ 
tifying localized peaks in particle density that are stable over 
multiple orbits. That is to say, an overdensity needs to have low 
radial drift and high longevity. For example, for Tf = 0.001 there 
are no strong overdensities, and for Tf = 3, the overdensities are 
short-lived and suffer strong radial drift. Therefore, our strategy 
is to integrate the particle density over a 25-orbit interval and use 
the Kolmogorov-Smirnov test to distinguish the resulting parti¬ 
cle distribution from a uniform distribution. We test several 25- 
orbit periods, spaced in steps of 10 orbits so that they overlap. 
The result of the KS test is a p-value, which measures the proba¬ 
bility that two observed data sets arose from the same underlying 
distribution. The method is discussed in Appendix A. We divide 
the results of the KS test into four categories: 


Clumping is unlikely: 
Clumping is somewhat likely: 
Clumping is likely: 

Clumping is very likely: 


p > 0.45 
0.25 <p< 0.45 
0.10 < p<0.25 

p<0.10 


Figure [^summarizes the result of the KS test for all the sim¬ 
ulations. The symbols in the figure mark the points where the 
test was performed. The figure is color -coded: For each point, 
the surrounding square has two colors (possibly the same color 
twice) that mark the extreme values of p. In this way, the figure 
indicates the range of values observed across the simulations. 
The figure shows that, within the streaming regime, there is a 
distinct region in the Z - Tf phase space where particle clumps 
are consistently likely to form. Of special interest is the fact 
that particles as small as Tf = 0.003 can form distinct particle 
clumps, at least some of the time, for solid concentrations around 
Z ~ 0.065. This value is high, but as we noted earlier, there is 
evidence that the ice line is associated with a particle concen¬ 
tration at least that high ( [Sirono|2011[ ), even before we consider 
turbulent eddies, disk winds or photoevaporation. 

Appendix B has convergence tests. Figure B shows the 
spacetime diagram of a 3D simulation with Tf = 0.03; Figure [B3] 




Fig. 2. Particle scale height for small particles (size Tf < 0.002, R < 1.5mm) over the course of the simulation. The left figure shows the entire 
simulation, including the sedimentation phase. After 250 orbits, has reached a steady state and the numerical experiment begins. The right 
figure shows only the gas removal phase. Gas is removed exponentially. In the right plot, the simulation time is replaced by the solid concentration 
Z = 2isoiid/^totai- For particles in this size scale, the scale height is determined primarily by A, and secondly by Z. 
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Fig. 3. Particle density Pp as a function of solid concentration (left) and height (right). The solid concentration is Z = Zsoiid/^^totai- The top figure 
shows Pp at the midplane (height z = 0) as a function of Z. The bottom figure shows pp at Z = 0.08 as a function of the z coordinate. The particle 
density is measured in terms of the initial gas density po- All the simulations have particle size Tf = 0.002. Two runs have a standard pressure 
support (A = 0.05) and two have a reduced pressure support (A = 0.025). Two runs have 128 x 128 grid resolution and two are 256 x 256. There 
is no evidence of the density cusp predicted b y|Sekiya|(|1998 1 and|Youdin & Shu| ( |2002| ). The A = 0.05 and A = 0.025 runs should have clear 
cusps after Z ~ 0.10 and Z ~ 0.02 respectively ( [Youdin & Shu 2002| ). These simulations may have an additional source of mixing that inhibits the 
formation of the Sekiya-Youdin-Shu cusp. The result is the same for smaller, Tf = 0.001 particles (not shown). 


shows the spacetime diagrams of higher resolution simulations 
with Tf = 0.003 and Tf = 0.010, and grid size 256 x 256. In all 
cases, the results are consistent with the spacetime diagrams in 
Figs. 0 an 

4. Constraints on planetesimal formation modeis 


[2009a| ). 


2. An alter native to high partic le concentration is low disk tur¬ 
bulence. Qrmel et al. ( 2008| ) have shown that, for a ~ 10“^, 
chondrule aggregates can reach ~ 4 - 5 mm. Figure 
shows how low a can be much more effective than high Z in 
triggering particle clumps. 


We consider the formation of dense particle clumps, as the pre¬ 
cursors of planetesimals. |Qrmel et~ar] ( |2008| ) have found that 
pm-sized dust sticks to chondrules to form a porous rim, which 
absorbs some of the kinetic energy from collisions. This allows 
the chondrules to form small aggregates, whose size depends 
primarily on the speed of the turbulent eddies in the disk, 


t^eddy 


; Va 


(23) 


where is the local sound speed, and a is a proportionality con 
stant that ni easures the strength of the disk turbulence. [Qrmel 


et al. 


( 2008| found that, for typical disk turbulence (a = 10“^), 
the chondrule aggregates have R ^ 1 mm in size. For low tur¬ 
bulence (a = 10“^), the chondrule aggregates have ~ 4 mm, 
with many aggregates as large as ~ 7 mm). 

For r > 10 AU, mm-sized aggregates have Tf > 0.01 and 
the streaming instability occurs easily. For r < 5 AU, Fig. [^im¬ 
poses some constraints on the plausible scenarios that can pro¬ 
duce dense particle clumps. These constraints are illustrated in 
Fig. In the figure, the “disk mass vs r” phase space is di¬ 
vided into regions where different disk models can form particle 
clumps. Broadly speaking, there are three avenues to form plan¬ 
etesimals in the inner solar system: 


3. Finally, since Tf is fundamentally a measure of the particle 
stopping time, as the disk clears and the gas density de¬ 
creases, Tf increases for all particles. Figure shows how 
the disk mass (M^isk) affects the location of the planetesimal 
forming region. The effect is most dramatic in the inner ~ 
1-2 AU, where particle clumping is not possible until M^isk 
has decreased. 

These scenarios are not exclusive. The formation of plan¬ 
etesimals may well rely on a combination of, say, disk evolution 
to decrease Mdisk and pressure bumps to increase Z. We wish 
to investigate which combination of (Z, a, Mdisk) are compatible 
with the formation of particle clumps in the disk. To do this, we 
generalize the Hayashi nebula to allow for different disk masses 


S = 1700 g cm 


-2 


/ r \-3/2 

/ Mdisk \ 

1 AU/ 

\ Mmmsn / 


(24) 


where Mmmsn is the disk mass for the Hayashi nebula. Equation 
[prelates Tf to E, the particle size R, and material density p,. Take 
a typical material density of p, = 3.6 g cm“^, and we obtain the 
minimum distance rpi where asteroids can form by the streaming 
instability, 


1. Figure shows that the streaming instability can be pushed 
to particles as small as Tf = 0.003 if the solids can reach 
sufficiently high mass fractions. There are various disk pro¬ 
cesses that may enhance Z sufficiently to induce clumping. 
They include the break-up or ice-dust agg regates (jSirono 
[201 Ij ), as well as large-scale pressure bumps ( [Johansen et al. 


Tpi = 4.33 AU 


Mdisk Tf.„in / R 

Mmmsn 0.003/ \ mm/ 


(25) 


where Tf^min is the minimum Stokes number needed to trigger the 
streaming instability, which is determined by the particle con¬ 
centration. In a disk where solids can concentrate to Z = 0.065, 
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Fig. 4. Spacetime diagrams showing the solid surface density Esom as a function of the radial coordinate ;r and simulation time. The right hand 
axis shows the mean solid concentration Z = (2isoiid)/<^totai)- The figure shows selected runs with Tf < 0.03. The top row has runs with standard 
pressure support (A = 0.05) and the bottom row has runs with a reduced pressure support (A = 0.025). The surface density is shown by a color 
scale. Some of the runs form visible filaments. For A = 0.05 there is a clear trend where lower Tf requires higher Z to form filaments. For 
A = 0.025 the behaviour is less predictable - while the particle concentration is higher, the low radial drift speed interferes with the streaming 
instability. 
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Fig. 5. Same as Fig.[^ but for larger particle sizes (Tf = 0.1 to 3). For Tf = 0.1 and A = 0.025 the simulation cannot resolve the minimum solid 
concentration (Z = Zsoiid/^^totai) that leads to clumping. All simulations begin at Z = 0.005 and that is sufficient to produce visible clumps for that 
Tf and A. Also, for Tf = 1 and A = 0.05, the filaments that form are short lived. 
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Solid concentration (Z) 


Fig. 6. Radial drift speed Vr of particles as a function of the particle 
size Tf and solid concentration Z = 2^soiid/^totai- Speed is measured in 
terms of the sound speed c^. All solids experience gas drag. Particles 
of size Tf ~ 1 experience very rapid drift, which decreases for larger 
particles. In addition, higher Z forces the gas component to orbit closer 
to the Keplerian speed, reducing |vr|. However, in all cases, the drift 
speed remains high. For reference, for |vr| = 0.005cs the particles cross 
the disk after a few thousand years. 


Fig. gives a minimum Stokes number of Tf^min = 0.003. Oth¬ 
erwise, take Tf^min = 0.006 and Z = 0.04. Equation (25) can be 
thought of as the location of the “planetesimal line”, or the min¬ 
imum distance where planetesimals can form. This line moves 
inward as the disk evolves. Figure shows the location of the 
planetesimal line for different choices of Z, a, and Mdisk- 


5. Coagulation for mm-sized particles 


Particles with R = 1 mm (rf = 0.003) may be a critical stage in 
the formation of planetesimals. These particles lie at the upper 
limits of coagulation, and the lower limits of the streaming insta¬ 
bility. The formation of planets may hinge on the ability of these 
particles to grow from Tf = 0.003, and move deeper into the 
streaming regime, where gravi tational collapse can occur (see 
Eq. [^and |Bai & Stone||2010b| ). In this section we explore the 
ability of mm-sized, Tf = 0.003 to grow faster than the disk ac¬ 
cretion timescale. 

Our simulations are too coarse to directly measure collision 
speeds between particles, but with some extrapolation, it is pos¬ 
sible to produce robust upper limits. To do this, we note that 
in a small region of the disk the particle velocities should be 
randomized and their relative speeds should follow a Maxwell- 
Boltzmann distribution, 


The Maxwell-Boltzmann distribution is parametrized by the 
characteristic speed a. Given a series of relative speeds 
the Maximum Likelihood Estimate for a is 


a = 


\ 3n 


1 " 
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(27) 




Radial coordinate x 


Fig. 7. Surface density of solids Esom as a function of the radial 
coordinate x, evaluated at two points in time 25 orbits apart (first blue, 
then green). The top plot is in the suspension regime (particle size Tf = 
0.001) and has solid concentration Z = (Zsoiid)/(^totai) = 0.095. In this 
regime, there is very little clumping even at high Z. The middle plot is 
the beginning of the streaming regime (Tf = 0.003) and Z = 0.095. In 
this regime, clumps form and are stable. The bottom plot is in the radial 
drift regime (Tf = 3) with Z = 0.05. In this regime particles form large 
over-densities that experience rapid radial drift. 


The strategy, then, is to extract a collection of particle pairs 
that are closer than some distance D and compute a from the 
relative speeds. This is shown in Fig. for several choices 
of D, and for three different simulations with Tf = 0.003. If 
we extrapolate a for D 0, we obtain a rough estimate of 
the Maxwellian distribution of collision speeds. If we fit a lin¬ 
ear relation, a = gq -r m D, the extrapolated value of a is 






















D. Carrera et al: How to form planetesimals from mm-sized chondrules and chondrule aggregates 



Fig. 8. Region of the particle size vs concentration phase space where the streaming instability is active. Particle size is measured in the stopping 
time Tf and the particle concentration is Z = Zsoiid/^^totai- The colors mark the probability that particle clumps can form, where red is “unlikely”, 
magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. See appendix A for the precise definitions and methods. When different 
simulations give different results, the two extreme values are shown. For example, if two simulations give “likely” and one gives “unlikely”, this 
is expressed as a square that is half red, half green. The symbols indicate the number of simulations available for the result. A triangle means 
three simulation, a cross means two, and a circle means one. Regions without symbols are extrapolations. The black lines indicate our (subjective) 
assessment of main regions where the behavior of particles is qualitatively different. 


ao = -0.0003 + 0.05 cm/s. Since the true ao cannot be nega¬ 
tive, we take the range 


0 cm s ^ < ao < 0.1 cm s 


(28) 


Figure [T^ shows the extrapolation of a for clq in this interval, 
and for the l-cr confidence region for the slope m. We are inter¬ 
ested in the probability that a collision of two mm-sized particles 
will result in sticking. For a material density of p, ~ 3.5 g cm“^, 
the critical speed is around Vcrit ~ 0.03 cm/s ( [Guttler et al.|2010| ). 
The probability that v < Vcnt is given by the cumulative distribu¬ 
tion function 


F(v) 


-f 


f(x) dx = erf 


<V2 


2 / -V 

^^ P |—2 


2oJ 


(29) 


Figure im shows the probability of a sticking collision F(vcrit) 
for the d values from Fig.fl^ As Z) ^ 0, the steady decrease in 
particle speeds leads to a significant increase in the fraction of 
collisions that result in sticking. 

Figure [T^ shows the Maxwellian distributions for ^ = 0.3 
cm/s (typical value for D = 10“^ H) and for d = 0.05 cm/s (mid¬ 
range value for ^o). The bulk of the collisions result in bouncing, 


and a small fraction result in sticking. For comparison, the figure 
also shows the root-mean-squared vertical speed Vnns for all the 
particles in the grid. The figure shows how Vnns greatly over¬ 
estimates the particle collision speeds. 

Finally, we use F’(vcrit) to determine the typical time needed 
before a 1 mm grain has a sticking collision. The key question 
is whether the particles can grow faster than the disk evolution 
timescale. The time needed for a sticking collision is given by 

4tick = [nilrfyup F(Vcrit)) *, (30) 

where is the particle number density, and v is the mean relative 
particle speed. In these simulations, ~ 4 x 10“^ m“^ at the 
midplane, and for a Maxwellian, v = 2a ^|2|7l. Figure [o shows 
the results for 4tick- As Z) ^ 0, we get 4tick ~ 10^ - 10^ years, 
which is a small fraction of the disk lifetime of ~ 10^ years. 
We repeated this experiment for several particle concentrations 
from Z = 0.01 to 0.1. Figure shows that a < 0.1 cm s“Ms a 
robust upper bound, largely independent of Z. We conclude that 
the particle growth timescale 4tick < 10^ years is also robust. 
Therefore, we conclude that Tf = 0.003 particles can double in 
mass within a small fraction of the disk lifetime. As the particles 
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Fig. 9. Location of the asteroid forming region in terms of disk mass 
Mdisk and semimajor axis r. This can also be thought of as a “disk age 
vs r” phase space. The phase space can be divided into three regions. 
In the red region (left), particle clumps cannot form for any disk model. 
In the yellow region (middle), particle clumps can only form if the disk 
turbulence is low, so that chondrule aggregate can reach ~ 4 - 5 mm 
in size (see main text). In the green region (right), particle clumps can 
form for ~ 1 mm particles and low turbulence is not required. Here, 
^MMSN is the mass of the minimum-mass solar nebula. The yellow 
and green regions can be further subdivided into an inner and outer 
sub-region. In the inner sub-region, particles have a stopping time of 
Tf > 0.003, and can only form clumps when the solid concentration is 
high - i.e. Z > Zsoiid/^^totai = 0.065. In the outer sub-region, the particle 
stopping time is Tf > 0.006 and clumps can form for Z > 0.04 (see 
Fig.[^. In all models, as the disk evolves, the disk mass drops and the 
planetesimal forming region moves closer to the star. 
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Fig. 10. Maximum likelihood estimate a of the characteristic speed 
between solid particles. The particle size is Tf = 0.003 and the solid 
concentration (Z = Zsoiid/^^totai) is 0.074 for all the simulations. The 
procedure is to select particle pairs with separation less than D and com¬ 
pute a from Eq. Extrapolating to D ^ 0 gives a rough estimate of 
the characteristic collision speed clq. The region in green corresponds to 
0 < ao < IcTa and the l-cr confidence interval for the slope. 


grow, the streaming instability becomes more effective and the 
density of the particle clumps increases. 
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Fig. 11. Probability that the relative speed between two solid particles 
lies below the sticking threshold Vcrit ~ 0.03 cm/s ( [Guttler et al.|2010| ). 
The particle size is Tf = 0.003 and Z = Zsoiid/^totai = 0.074 in all the 
simulations. Given a set of particle pairs with separations less than D, 
one can fit a Maxwellian distribution (Figs. and determine the 
fraction of speeds below Vcrit- Extrapolating to D ^ 0 gives an estimate 
of the fraction of particle collisions that result in sticking. The region in 
green corresponds to the l-cr confidence region marked in Fig.fTOj 


6. Conclusion 

The formation of planetesimals remains a difficult problem. The 
traditional bottom-up coagulation model faces important chal¬ 
lenges from theory and observation. There is growing evidence 
that planetesimals must form top-down , from the gravitationa l 
collapse of dense clumps of solids (e.g. jMorbidelli et al.|2009| ). 
In this paper we present the results from over 100 hydrodynamic 
simulations where we modeled the dynamics of solid particles 
inside a protoplanetary disk. We model particle sizes from sub¬ 
millimeter-sized chondrules to meter-sized rocks, and evaluate 
the effect of particle concentration and gas pressure gradient. 
Our key results can be summarized as: 


1. We present the particle properties as a size vs concentration 
phase space, and we map the region where particle clumps 
can be expected to form (Fig. [^. We measure the particle 
size as Tf = where tf is the friction time and Ok is 

the Keplerian frequency. We measure particle concentration 
as Z = ^soiids/^totah where E denotes surface density. We 
find that particles as small as Tf = 0.003 can participate 
in the streaming instability and form dense clumps at 
Z ~ 0.065. Although this solid concentration is quite high, 
it is smaller than what can be produced by the break-up of 
ice-dust aggregates near the snow line ( |Sirono|2011| ). More 
generally, high Z values may be achieved by the radial drift 
of solids from the outer disk into the inner disk ( jYoudin &| 
Shu||2002|), and by large scale pressure bumps or vortices 


( [Johansen et al.|2009a| ). 


2. We find constraints on asteroid formation models and 
present them in Fig. We map the region in the Mdisk 
vs r phase space that is consistent with the formation of 
dense particle clumps. The location of the asteroid-forming 
region is primarily a function of disk mass and turbulent 
viscosity, as low a leads to larger chondrule aggregates 
Qrmel et~ar] ( [2008[ ). After that, the most important factor 
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Fig. 12. Maxwellian distribution of relative speeds for mm-sized particles (Tf = 0.003, pressure support A = 0.05 and solid concentration 
^ - ^soiid/^totai = is 0.074). For particles separated by distance d < 10“^// {H is the disk scale height) the Maxwellian distribution has characteristic 
speed a ~ 0.3 cm/s (red). Extrapolating as d 0 we obtain a con servative estimate of do = 0.05 cm/s (blue). For reference, we include the 
sticking, bouncing and fragmenting regions from |Guttler et al.| ( |2010| ), and the root-mean-squared speed Vnns for all the particles in the simulation. 
Note that the shape of the curve is distorted by the log scale; in particular, the great majority of collisions are in the bouncing regime. Also note 
that Vnns is a poor estimate for the the collision speeds. The value d is obtained by the maximum likelihood method. 
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Fig. 13. Time between sticking collisions. The figure shows the 
three 128^ runs from Fig.|ll|plus two higher-resolution runs. In all runs 


the particle size is Tf = 0.003 and Z = 2^soiid/^totai = 0.074. We fit a 
Maxwellian distribution to to all particle pairs with separation less than 
D. Extrapolating the fit toward small D we estimate the time between 
sticking collisions. The points on the figure appear to follow a power 
law 4tick ^ The r egio n in green corresponds to the l-cr confidence 


region marked in Fig. 10 The figure shows that 4tick <10^ years is a 


robust upper bound on 4tick- 


Fig. 14. Characteristic collision speed do between solid particles at 
different points in our simulations, with l-cr error bars, for particle size 
Tf = 0.003. The solid concentration (Z = Zsoiid/^^totai) has little effect on 
do. The value ao < 0.1 cm/s (red dashed line) is a robust upper bound 
on do. This corresponds to coagulation timescale of 4tick <10^ years, 
which is significantly less than the disk lifetime (10^-10^ years). Note 
that the fitting method can give negative values because it is a simple 
linear extrapolation (see main text). 


is the disk’s ability to increase the solid concentration locally. 
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3. We estimate the probability of sticking collisions between 
mm-sized particles (rf = 0.003). We find a robust upper 
limit on the timescale for these particles to stick and grow to 
larger sizes (~ 10^ years) that is significantly smaller than the 
lifetime of the disk (~ 10^ years). A high concentration of 
particles is not required for this result. As the particles grow, 
the streaming instability becomes more effective, leading to 
the formation of planetesimals. 

Altogether we find that particle concentration by the stream¬ 
ing instability provides a viable path to forming asteroids di¬ 
rectly from mm-sized chondrules, particularly if weak turbu¬ 
lence facilitates the growth of chondrule aggregates with sizes 
of a few millimeters. 
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Appendix A: A measure of particle clumping 

Particle clumping is difficult to measure objectively. Previous 
authors have measured the peak particle density, but this is a 
very local measure that contains limited information about the 
particle behavior. In particular when small particles (rf < 0.003) 
form filaments that are clearly visible in a spacetime diagram 
(Fig. 15, no single grid cell achieves a very high density. At the 
same time, large particles (rf > 3) have may have a few cells 
with very high density, but the density peaks are unstable and 
short lived. 

Figure [ATT] illustrates our approach. For each time step, we 
take the particle surface density Ssoiid (top row). Particle over¬ 
densities are clearly visible on this plot. We then average Ssoiid 
over 25 orbits ((2soiid)t, second row). This serves to to smooth 
out overdensities that are unstable, short-lived, or have too much 
radial drift. We chose 25 orbits because it is half of the 50-orbit 
timescale for Z to increase. We considered using the radial speed 
to separate clumps that are drifting into the star. The appeal of 
radial speed is that it provides an instantaneous measure. How¬ 
ever, radial speed does not capture the longevity of a particle 
clump. For example, as in a traffic jam, the pattern speed of the 
clump may be lower than the particle drift speed. At the same 
time, particles with low radial speed may produce short-lived 
clumps. Therefore, it is necessary to compare the location of a 
clump across time. Averaging Ssoiid over time is a simple way to 
achieve this goal. 

Finally, we sort the grid cells of <2soiid)t in order from highest 
to lowest density, and compute the cumulative distribution F of 
the sorted grid cells. If the particle distribution is uniform, then 
F forms a straight diagonal line G from (0,0) to (1,1). We use 
the difference F - G as a large-scale measure of clumping. The 
Kolmogorov-Smirnov ( |Smirnov|1948[ Wall & Jenkins |20()3] ) test 
gives a standard way to compare F and G. It computes a p- 
value that measures probability that the data set that produced F 
follows the distribution G 


Qiz) = 2j^{-iy-Uxp{-2z^f) (A.1) 

p = Q{D^\ (A.2) 

where D is the maximum distance between F and G, and n is 
the number of independent observations. We caution the reader 
against interpreting these p-values as strict probabilities; we only 
use as a convenient metric to measure clumping. With that in 
mind, we divide the p-values into four broad categories: 
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Clumping is unlikely: p > 0.45 

Clumping is somewhat likely: 0.25 < p < 0.45 

Clumping is likely: 0.10 < p < 0.25 

Clumping is very likely: p < 0.10 

Figure shows the result of this analysis across our entire 
range of simulations. The figure is broadly consistent with the 
spacetime diagrams in Figs. an d 0 

Appendix A. 1: Alternate measures of particle clumping 

To confirm our results, we devised an alternate method to mea¬ 
sure particle clumping, based as much as possible on very differ¬ 
ent principles. Our second method is to perform a chi-squared 
test on {pp}zu without sorting the grid cells. We use the reduced 
chi-squared statistic. 


2 ^ ^ ^ i V ~ 
^ ~ 0-2 


(A.3) 


where k = n - 1 is the degrees of freedom, Oj is the obser¬ 
vation (value of {pp)zt on grid cell 7 ), Ej = (Pp) is the ex¬ 
pected density for a uniform distribution, and is the expected 
variance between different measurements. > 1 , that indi¬ 

cates that the difference between observations and the model is 
greater than can be explained by the variance cr^. This can be 
because the model is incorrect, or because the true variance be¬ 
tween measurements is greater than cr^. We divide into four 
intervals: 

Clumping is unlikely: < 0-7 

Clumping is somewhat likely: 

Clumping is likely: 1.0 < AfJed 

Clumping is very likely: x^q^ >1.3 


That leaves as the only unspecified parameter. When 
the simulations start, is almost zero, giving unrealistically 
high;^^^^ values. As sedimentation proceeds, the value of in¬ 
creases. We opted for a that is typical for the last 25 orbits 
of the sedimentation phase. Becauseis sensitive to cr^, there 
is necessarily some amount of subjectivity in cr^. In the end, 
we chose so that the Tf = 0.003 simulations wold give the 
same result as in Fig. Therefore, the true test is whether the 
xl^^ method remains consistent with the KS method across all 
the other p artic le sizes. 

Figure |A.2| shows the final results for the x^^^ method. The 
similarity between this figure and Fig. is strilang. Consider¬ 
ing how different the two techniques are, the agreement between 
Figs. [ 8 ] and I A. 2| indicates that our core results are robust. 


where u is the gas speed, g = ffz is the vertical gravitational 
acceleration, and p is the effective fluid density (treating gas and 
solids as a single fluid). |Bai & Stone] pOlO'al have shown that 
the streaming instability is able to keep the Richardson number 
above the critical value needed for the KHI, so that the KHI is 
absent in this problem. To confirm this, we performed two 3D 
simulations with Tf = 0.03 and A = 0.05. The experiment per¬ 
formed in our 3D simulation is slightly different from that of our 
2D simulations. Instead of increasing Z gradually, we only test 
the behavior of the system at three discrete values: Z = 0.01, 
Z = 0.02 and Z = 0.03. The reason for this is that 3D simu¬ 
lations are very computationally expensive, especially for small 
particles, so that a long-term simulation is prohibitive. Our 3D 
runs have three phases: 


- Phase /: Each run starts with Z = 0.01, and Z is held constant 
for 30 orbits. 

- Phase IP. The value of Z is increased quickly. In run “3D.Z2” 
we increase Z to Z = 0.02, and in run “3D.Z3” we increase 
ZtoZ = 0.03. 

- Phase IIP. The value of Z is once again held constant for the 
remainder of the run. 


Figure [BT] shows the spacetime diagrams for both runs. The 
figure shows that both runs produce visible particle clumps for 
Z > 0.02 and no clumps for Z = 0.01. It is also clear that 
clumps form more readily with Z = 0.03. These results are con¬ 
sistent with our 2D simulations (Fig. |^. For particles smaller 
than Tf = 0.03, the small integration steps mak e 3D simulations 
prohibitive. However, as explained in section |2.6[ the long re¬ 
sponse time of very small particles (Tf < 0.0064) means that our 
simulations will overestimate Z; so our results can be considered 
robust upper bounds for this size range. 

Figure [R2] gives a side view {x-z axis) and a top view {x-y 
axis) of the 3D.Z2 run dXt = 150 orbits. The clumps and their 
filamentary structure are very prominent. An important feature 
of the 3D runs is that they form visible clumps more easily than 
the 2D runs. Figure [B3] shows the formation of clumps in the 
two 3D runs and one of the 2D runs. This figure shows that the 
3D runs are consistent with the 2D run. In particular, the clumps 
are readily visible at Z = 0.02 and the 3D runs show no evidence 
of additional stirring from the KHI. 

Figure IB. 41 shows the maximum particle density (Pp,max) for 
the two 3D simulations and one of the 2D simulations. One 
salient feature of the figure is that in the early phase of the sim¬ 
ulations, before clumping occurs, the 3D runs consistently show 
a peak density 2-3 times greater than the 2D runs. This occurs 
because particles can also accumulate along the azimuthal direc¬ 
tion. The figure also shows Pp,max for the 3D runs after averaging 
along the azimuthal direction (purple lines). This averaging re¬ 
moves the apparent discrepancy between the 2D and 3D runs. In 
the later stages the 2D runs reach higher peak densities because 
they have very high Z values. 


Appendix B: Convergence tests 

Appendix B. 1: Three-dimensional box 

Our simulations are axisymmetric, they neglect the effect of the 
Kelvin-Helmholtz Instability (KHI). Whether the KHI is present 
is dictated by the Richardson number ( [Chandrasekhar| 19^ , 


Appendix B.2: Resolution 

We performed two simulations at 256 x 256 grid resolution with 
Tf = 0.003 particles. Figure [B3] shows the spacetime diagrams 
for both sets of simulations. Since the behavior of the 256 x 
256 runs is almost identical to the corresponding 128 runs, we 
conclude that, in terms of resolution, our simulations are fully 
converged. 
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Fig. A.l. Particle density distribution for three particle sizes. Particle size is measured in the stopping time Tf, where Tf = 0.001, 0.003 and 
3 correspond to the suspension regime, streaming regime and radial drift regime respectively. Row (a) is the Surface density of solids Esom as a 
function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). Row (b) averages Esom over the 25-orbit 
interval (Esoiid)t- Row (c) shows (Esoiid)t with the grid cells sorted from highest to lowest density. A steep curve indicates stable particle clumps. 
Row (d) shows the cumulative distribution of (Esoiid)t in row (c) (solid red curve) along with the cumulative distribution of a perfectly uniform 
particle distribution (dashed blue line). The distance between these two curves serves as a measure of particle clumping. In these plots the particle 
concentrations Z = (Esoiid)/(^totai) are (left to right) 0.095, 0.095 and 0.05. 
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Fig. A.2. As in Fig.j^ but using the;^^^^^ method to estimate the likelihood of clumping instead of the KS-derived method. The figure marks the 
region of the particle size vs concentration phase space where the streaming instability is active. Particle size is measured in the stopping time Tf 
and the particle concentration is Z = Zsoiid/^totai- The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is 
“somewhat likely”, blue is “likely”, and green is “very likely”. When different simulations give different results, the two extreme values are shown. 
The symbols indicate the number of simulations available. The circle, cross and triangle indicate, respectively, one, two and three simulations. 
Regions without symbols are extrapolations. To facilitate the comparison with Fig. the boundary marked by the black lines has been copied 
from Fig. [^without modification. 
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Fig. B.l. Spacetime diagram of two 3D runs with particle size Tf = 0.03 and A = 0.05. The color indicates the column density Zsoiid/(2^soiid) 
using the same color bar as Figs.an d0 Both runs begin with Z = 0.01 and have Z increased only for a short interval. In one run (left), Z grows 
to 0.02 and in the other (right) Z grows to 0.03. Clumps are visible for Z > 0.02, consistent with the 2D runs. 
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Fig. B.2. A snapshot of run 3D.Z2 taken at ? = 150 orbits (Z = 0.02). Image (a) is the view from the side (x-z plane). Image (b) is the view 
from above (x-y plane). The color scale marks the column solid density normalized to the initial gas density (Zsoiid/^gas,o)- The two images have 
a different color scale, since the column density in the azimuthal direction is much higher. Note that the particle structure is nearly axisymmetric. 
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Fig. B.3. A sequence of vertical snapshots (x-z plane) showing the formation of distinct particle clumps in one of our 2D simulations (left) and 
our two 3D simulations (middle, right). For each run, the images are taken 10 orbits apart. The time and Z values are shown in each plot. The 
times were chosen to show the formation of the first distinct clumps, and to give a simil ar pe ak density on the final image (pp^max = 1-8,1.5,2.1 
left to right). The color indicates the column density using the same color scale as in Fig. |B.2[ a). The other 2D runs give similar results to the one 
shown. The 2D run is consistent with the 3D runs. In particular, the 3D runs do not show any evidence of additional turbulence stirring compared 
to the 2D runs. 
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Integration time after sedimentation (t/27rfl ^) 

Fig. B.4. Maximum particle density (Pp,max) for one of the 2D simulations and the two 3D simulations. The top plot compares run 3D.Z2 (blue) 
and a 2D run (red). The bottom plot compares run 3D.Z3 (blue) and the same 2D run (red). In each plot the purple line shows Pp,max for the 3D run 
after averaging along the azimuthal direction. The vertical dashed lines mark the places where the 3D runs transition from Z = 0.01 to Z = 0.019 
(for 3D.Z2) or Z = 0.029 (for 3D.Z3). Similarly, the vertical solid lines mark the progressive increase of Z during the course of the 2D mns. 
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Fig. B.5. Spacetime diagrams showing the solid surface density Esom as a function of the radial coordinate x and simulation time for Tf = 0.003 
and Tf = 0.010 particles. The surface density is shown by color, using the same color scale as in Figs.|^an The figure shows both 128 x 128 
simulations and the corresponding 256 X 256 simulations. In all simulations A = 0.05. The higher-resolution runs produce results consistent with 
the 128 X 128 runs. 
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